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FOREWARD 


This  technical  report  was  prepared  by  Drs.  William  Baumann,  Terry  Herdman,  Harold 
Stalford,  and  Mr.  Frederick  Garrett,  Jr.  This  work  was  performed  under  the  sponsorship 
of  the  United  States  Air  Force  (USAF),  and  was  administered  under  the  direction  of  the 
Air  Force  Flight  Dynamics  Laboratory,  Air  Force  Systems  Command,  Wright-Patterson  .Air 
Force  Base,  Ohio.  Mr.  Charles  F.  Suchomel  (AFWAL/FIGC)  was  the  Project  Engineer 
for  the  Air  Force.  This  report  represents  a  part  of  the  task  to  identify  a  technique  that 
describes  nonlinear  aircraft  dynamic  properties  potentially  useful  as  flying  quality  indicators 
in  large-amplitude  and/or  high-acceleration  motions.  This  work  was  accomplished  over  the 
period  16  July  1986  through  15  July  1987. 


1.  Introduction 

High  performance  aircraft  have  mission  requirements  for  operating  in  high  angle-of-attack 
and  sideslip  conditions  and  in  various  maneuvers  of  unsteady,  large  amplitudes  governed  by 
multiple  input,  multiple-output  nonlinear  dynamics.  The  large  amplitudes  coupled  with 
high  degree  nonlinear  dynamics  lead  to  requirements  for  nonlinear  flying  qualities.  These 
qualities  require  analyses  of  nonlinear  dynamic  stability,  nonlinear  coiitroh  and  nonlinear 
response  behavior  of  such  aircraft.  Techniques  are  needed  that  (1)  are  computationally  fea¬ 
sible,  (2)  retain  the  nonlinearities  of  the  aircraft  and  (3)  provide  physical  and  mathematical 
understanding  of  the  important  nonlinear  flying  qualities. 

One  methodology  which  has  recently  exhibited  promising  evidence  of  accomplishing  the 
above  is  the  Volterra  Series.  The  purpose  of  this  project  is  to  obtain  Volterra  series  rep¬ 
resentations  of  nonlinear  systems  typical  of  high  performance  aircraft  in  large  amplitude 
maneuvers. 

The  development  of  the  Volterra  representation  for  the  six  DOF  (degree  of  freedom) 
rigid  body-equations  of  motion  is  given  in  [1].  Therein,  the  nonlinear  equations  of  motion 
are  formulated  using  the  four-parameter  (quaternions)  method  (see  [2,3,4]).  The  first  three 
terms  of  the  Volterra  Series  representation  are  constructed  in  [5]  for  the  longitudinal  motion 
of  the  F-8  crusader  aircraft  using  the  nonlinear  model  developed  by  Gerrard  and  Jordan 
[6].  Suchomcl  [5]  has  numerically  evaluated  the  three- term  Volterra  Series  representation 
of  an  extended  version  of  Garrard  and  Jordan’s  model  which  incorporated  some  additional 
aerodynamic  drag  and  thrust  terms.  Suchomel’s  work  compares  the  predictive  performance 
of  the  three-term  Volterra  Series  representation  against  the  Runge-Kutta  integration  of  the 
full  nonlinear  model.  His  results  demonstrate  that  the  approximating  three- term  Volterra 
Series  gives  an  accurate  prediction  of  the  output  response  of  the  given  nonlinear  system. 

The  present  work  continues  the  application  of  the  Volterra  Series  representation  [1]  to 
nonlinear  aircraft  models  of  high  q  flight.  The  nonlinear  models  of  such  flight  have  natural 
low-order  submodels  which  hold  over  separate  flight  regimes  (e.g.,  pre-stall,  stall  and  post- 
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stall).  For  example,  ia  the  pre-stall  regime  we  observe  an  almost  linear  relationship  between 
the  state  variables  and  the  aerodynamic  force  and  moment  coefficients.  In  the  stall  and  post¬ 
stall  regimes  the  relationship  becomes  a  combination  of  linear,  bilinear,  quadratic,  cubic,  etc. 
It  is  not  only  necessary  but  unnatural  to  search  for  a  single  high-order,  nonlinear  relationship 
that  governs  flight  behavior  in  the  total  envelope  spanning  pre-stall,  stall  and  post-stall  flight. 
It  is,  however,  feasible  to  provide  simple  low-order  dynamics  that  govern  each  separate  flight 
regime  and  then  combine  them  to  form  the  total  nonlinear  model.  The  total  then  consists  of 
a  set  of  simple  low-order  equations  that  govern  flight  in  each  regime  (or  subspace)  and  that 
agree  with  the  equations  of  the  adjoining  regime  at  their  common  boundaries.  In  this  work 
we  consider  such  subspace  models.  They  come  about  naturally  as  described  above. 

In  Section  3.1  we  consider  a  nonlinear  wind  tunnel  model  for  high  alpha  (a)  longitudinal 
flight  involving  the  limit  cycle.  We  show  how  the  complex  nonlinear  aerodynamic  model  has 
a  natural  representation  in  terms  of  linear  and  quadratic  subspace  models  which  have  fairly 
simple  Volterra  Series  representations.  Simulation  analysis  is  conducted  on  both  the  original 
nonlinear  windtunnel  model  and  the  Volterra  subspace  approximation.  The  accuracy  of  the 
approximation  is  investigated. 

Next,  in  Section  3.2,  we  consider  a  nonlinear  windtunnel  model  of  wing  rock  which 
is  generated  by  unsteady  aerodynamics  effects  [23-33].  The  complex  model  is  a  simple 
composition  of  two  bilinear  subspace  models.  Again,  simulation  analysis  is  conducted  and 
the  accuracy  of  the  approximation  is  investigated. 

Section  4  considers  methods  for  analyzing  the  nonlinear  responses  discussed  in  Section  3. 
In  4.1  an  approach  based  on  Volterra  models  is  discussed.  The  application  of  classic  analysis 
techniques  to  the  sanie  models  of  Section  3  is  explored  in  Section  4.2. 

Our  conclusions  and  suggestions  for  future  research  are  contained  in  Section  5. 
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2.  Volterra  Series  Representation  For  Nonlinear  Systems 


The  first  approach  usually  taken  when  confronted  with  the  task  of  analyzing  a  nonlinear 
system  of  ordinary  differential  equations  is  to  linearize.  However,  some  current  flight  dynam¬ 
ics  problems  of  interest  to  the  Air  Force  possess  significant  nonlinearities  in  the  complete 
mathematical  model  that  cannot  be  ignored.  To  predict  the  behavioral  characteristics  of 
the  system,  the  complete  nonlinear  model,  or  at  least  an  approximation  that  contains  the 
significant  nonlincMrities,  must  be  analyzed. 

The  well-developed  theory  of  differential  equations  has  provided  general  existence  and 
uniqueness  results  for  nonlinear  systems.  However,  these  results,  unlike  the  corresponding 
results  for  linear  systems,  do  not  provide  a  closed  form  representation  for  the  solution. 
The  mathematical  tools  for  the  analysis  of  nonlinear  systems  are  limited,  and  the  known 
techniques  are  not  well  developed  when  compared  to  the  results  for  linear  systems. 

One  technique  for  the  study  of  nonlinear  systems  that  has  shown  promise  is  the  Volterra 
series  approach  [1,5,7-19].  This  approach  provides  a  series  representation  for  the  input-output 
behavior  of  the  nonlinear  system.  In  particular,  the  Volterra  series  gives  a  mathematical 
representation  of  the  solution  in  the  form  of  an  expanding  infinite  series  of  integrals  which 
encompasses  the  nonlinearities  of  the  system.  This  representation  can  be  viewed  as  a  series 
in  which  each  term  is  the  solution  of  a  linear  equation  where  the  nonlinearities  of  the  system 
appear  as  forcing  functions.  For  a  nonlinear  differential  equation  of  the  form 

X  =  fix,u) 


the  Volterra  series  representation  of  the  solution  is  an  infinite  series  of  integrals 

x{t)  =  ho{t)  +  /  /ji(f  - 
Jo 

+  /  /  h2{t  -  —  a2)u{ai)u{a2)d(Tida2 

Jo  Jo 

rt  ft  ft 


(2. 


Iff 


h:i{t  —  ai,  <  —  (72,  t  —  a3)u{ai)u{a2)u{(T3)da\da2da'i 
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where  ho{t)  is  the  zero-input  solution  and  the  /i,, ?  =  1,2,...,  are  called  the  ith  order  Volterra 
kernels.  For  a  complete  discussion  of  the  Volterra  series  representation  for  ordinary  differ¬ 
ential  equations  see  [20,21].  Also,  see  Appendix  A  for  additional  comments  on  the  \blterra 
series. 

The  Volterra  series  approach  provides  an  approximation  for  the  solution  by  truncating 
the  series.  Such  a  finite  representation  gives  the  response  of  the  system  in  terms  of  the  input 
signals,  initial  conditions,  and  the  stability  derivatives  of  the  system.  This  approximation 
provides  not  only  a  methodology  for  predicting  the  behavioral  responses  of  the  nonlinear  sys¬ 
tem  but  also  a  methodology  for  the  study  of  how  various  parameters  affect  the  controllability 
and  stability  of  the  system. 

The  mathematical  theory  for  the  Volterra  representation  is  by  no  means  complete.  .M- 
though  some  convergence  results  for  the  series  representation  exist,  they  are  extremely  con¬ 
servative  and  have  not  proven  useful  in  engineering  applications.  In  most  applications  the 
.■ligiiincaiiL  ^’iiciiacteristics  of  the  system  are  contained  in  the  second  or  third  terms;  therefore 
convergence  is  not  a  major  issue. 

Our  technique  for  approximating  the  solution  of  nonlinear  models  of  high  a  flight  is 
based  on  a  subspace  concept.  The  subspaces  arise  naturally  in  the  aerodynamic  forces  .and 
moments  in  high  performance  flight.  Our  approach  is  outlined  as  follows: 

1.  partition  the  total  envelope  into  natural  subspaces  in  which  the  nonlinearities  can  be 
accurately  described  by  low-order  multivariable  polynomials  (preferably  third  order  or 
less), 

2.  pick  an  equilibrium  point  for  each  subspace,  (we  have  no  definite  procedure  for  selecting 
the  equilibrium  point  -  at  present  these  points  are  selected  in  an  ad  hoc  manner), 

.3.  approximate  the  solution  of  the  differential  equation  for  each  subspace  by  a  truncated 
Volterra  series  (three  terms  or  less). 

For  this  approach  the  differential  equation  form  is  employed  for  the  terms  of  the  Volterra 


4 


series  instead  of  the  integral  from  (see  Appendix  A).  Consider  inputs  of  the  form 

u  =  uo  + 

where  uq.  U]  denote  an  equilibrium  and  a  perturbation  input.  The  solution  to  the  nonlinear 
differential  equation  has  the  form,  see  [20,  Chapter  three], 

X  —  Tq  kx\  -j-  k^X2  +  k^x^  -j_  . . . 

where  x  is  the  state,  k  is  a  constant,  Xo  is  the  initial  equilibrium  state  and  x,  is  the  contribu¬ 
tion  to  the  state  from  the  ith  term  in  the  Volterra  series.  Substituting  this  representation  for 
the  solution  x  into  the  governing  equation,  expanding  the  right  hand  side  in  a  Taylor  series 
with  respect  to  k,  and  equating  coefficients  of  equal  powers  of  k  gives  ri.se  to  the  following 
equations  for  the  first  three  terms  of  the  Volterra  series  representation 

x\  =  Axi  +  6ui 

i-2  =  Ax2 +  5ri(xi),  X2(0)  =  0 

X3  =  Ax3  + g2{xi,X2),  13(0)  =  0 

where  .4  is  a  n  x  n  matrix  (the  state  is  n  x  1),  b  is  a  n  x  1  vectui,  ^nd  the  y,,  f  =1,2  are 
in  general  n  x  1  nonlinear  functions  of  their  variables.  All  of  these  quantities  depend  on  the 
equilibrium  point  about  which  the  Volterra  series  is  expanded.  Upon  entering  a  subspace  the 
state  equations  are  reinitialized  by  incorporating  the  current  state  into  the  initial  condition 
of  the  linear  part  xy.  In  the  event  that  the  polynomial  of  step  1  above  is  first  order,  the 
corresponding,  gi,i  =  1,2  are  zero  and  the  Volterra  series  is  simply  a  linear  system. 

In  .Appendix  G  we  describe  a  program  that  computes  the  differential  equations  for  each 
Volterra  term.  This  program  was  used  for  all  the  Volterra  equation  calculations  discussed  in 
this  report. 
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3.  Simulation  Results 

The  following  examples  will  serve  to  illustrate  the  application  of  the  Volterra  series 
representation  to  subspace  models  of  nonlinear  systems. 

3.1  Longitudinal  Limit  Cycle:  Example  1 

For  the  purpose  of  illustrating  our  approach,  we  consider  a  simplified  nonlinear  model  of 
the  longitudinal  limit  cycle  at  high  angle  of  attack,  a.  Simplified  high  q  flight  is  governed 
roughly  by  the  following  two  differential  equations  developed  in  Appendix  B 

d  =  <7  +  9.168C(q)  -  1.8336((5,  +  T)  +  7.361904  (3.1 ) 

q  =  5.73(C„,q  +  Cmje)  +  2.865  (3.2) 

where  a  is  the  angle-of-attack  in  degrees,  q  is  the  pitch  rate  in  degrees  per  second  and  8^ 
is  the  elevator  control  in  degrees.  To  keep  the  presentation  simple,  we  prescribe  6’^.,  =  —1 
and  =  —1.5.  The  nonlinear  plunging  force  coefficient  C^{a)  is  represented  by  the 

following  equations  (3.3)  through  (3.6).  It  h^ls  the  appearance  of  an  inverted  high  o  lift 
curve,  Figure  3-1.  This  model  In  discrete  data  points  was  taken  from  measured  wind-tunnel 
values  of  the  T-2C  airplane,  Stafford  [22).  The  coefficients  in  equations  (3.3)  through  (3.6) 
were  numerically  calculated  to  provide  a  fit  for  the  discrete  data  points  for  the  corresponding 
Q  intervals. 

The  slope  of  the  Cj  curve  is  approximately  linear  up  to  stall  which  occurs  around  a  = 

.  This  portion  of  the  Cj  curve  is  represented  accurately  by  the  linear  function 

Cffo)  =  -G.07378494Q,  a  <  14.36T  (3.3) 

In  the  stall  region  between  14.36^'  and  15.6'’  the  curve  is  quadratic  in  nature  and  is 
represented  accurately  by  the  quadratic  equation 

C'ffo)  =  0.097220-*  -  2.86.530  +  20.03846,  14. .36°  <  o  <  15.6°.  (3.4) 

In  the  stall/post-stall  region  between  15.6'’  ainl  19.6°  the  C,  curve  revf'rses  its  curvature'  arul 
is  re[)resented  accurately  liy  the  quadratic  e<|uation 

(",(o)  =  -0.01971 0-*  +  0.74391  o  -  7.807.53,  15.6°  <  o  <  19.6°  (3.5) 


fi 


o  (DEGREE) 


REGION  I;  o  <  14.36°,  ~  LINEAR  FUNCTION 

REGION  II:  14.36°  <  q  <  15.6°,  C,(q)  ~  QUADRATIC  FUNCTION 
REGION  III:  15.6°  <  q  <  19.6°,  C,(q)  ~  QUADRATIC  FUNCTION 
REGION  IV:  19.6°  <  a  <  28°,  r,(a)  ~  LINEAR  FUNCTION 
FIGURE  3-1.  NONLINEAR  PLUNGING  FORCE  COEFFICIENT  C,(ft) 
FOR  T-2C  AIRPLANE  =  -IS° 


In  the  post-stall  region  between  19.6°  and  28°,  the  Cz  curve  exhibits  linear  behavior 

C,(o)  =  -0.01667Q  -  0.47333,  19.6°  <  a  <  28°.  (3.6) 

The  above  division  of  the  a  interval  into  four  a  subspaces  leads  to  an  accurate  representation 
of  the  nonlinear  function  Cz{oi)  by  four  low-order  equations. 

In  each  of  the  subspace  regions  described  by  Equations  (3.3)  through  (3.6)  an  equilibrium 
point  is  chosen  and  the  solution  of  the  differential  equation  is  approximated  by  the  first  three 
terms  of  its  Volterra  series.  For  simulation  and  analysis  purposes,  we  use  the  differential 
equation  form  for  the  N'olterra  series  terms.  Following  the  procedure  outlined  in  Section  2, 
with  state  .r  =  and  input  u  =  8^.  (elevator  angle),  we  consider  inputs  of  the  form 

8,  =  8,,  -I-  k8,,  (3.7) 

where  8e^  is  the  equihbrium  input  and  k  is  a  real  parameter  needed  for  the  perturbation  anal¬ 
ysis  below  and  does  not  have  physical  significance.  The  solution  to  the  nonlinear  differential 
equations  (3.1)  and  (3.2)  has  the  form 

O'  =  Xo  +  kxx  +  k’^X2  +  k^X:i  +  •  •  •  (3.8) 

where  Jq  =  [ao,<7o]^  is  the  equilibrium  state  and  x,  =  is  the  contribution  to  the  state 

from  the  ith  term  in  the  Volterra  series.  Substituting  (3.7)  and  (3.8)  into  (3.1)  and  (3.2), 
expanding  in  a  Taylor  series  with  respect  to  k^  and  equating  coefficients  of  equal  powers  of 
k  results  in  differential  equations  for  the  first  three  Volterra  series  terms  of  the  form  ^ 


ii  =  Axi  +  b8f.^ 


i2  =  8Kx2^- gx{xx),  X2(0)  =  0  (3.9) 

X3  =  ylx3  -I-  g2{xx ,  X2),  13(0)  =  0. 


In  all  four  subspaces  regions  the  column  matrix  b  is  given  by  (sec 


.Appendix  C  for  Region  11) 


-1.8336 

-8.695 


(3.10) 
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and  the  matrix  A,  the  Jacobian  at  (oq, 9o, ^eo)»  has  the  form 

-5!73  0 

The  entry  /In  depends  in  this  example  on  the  value  oo-  However,  in  subspace  regions  I 
and  IV,  as  described  by  equations  (3.3)  and  (3.6),  the  approximation  (3.9)  is  linear  with 
(5rj(xi)  =  g2{xi,X2)  =  0  and  An  is  independent  of  ao  with  values  —0.67646  and  —0.152831, 
respectively. 

In  the  subspace  regions  II  and  III,  as  described  by  Equations  (3.4)  and  (3.5),  the  functions 
Qi  and  g2  have  the  form 

(3.12) 

g2{Xl,l2)=  ^21^1  <^2  (3^3) 

where  the  constants  g^  and  5^21  depend  on  the  particular  region.  We  picked  the  equilib¬ 
rium  point  for  each  region  by  first  assigning  a  value  for  6^^  and  then  calculating  ao,7o  from 
equations  (3.1)  and  (3.2).  In  region  II,  the  Volterra  series  has  been  expanded  about  th  *  equi¬ 
librium  point  (ao5  9oi<5eo)  =  (14.36'’, —1.7552'’/sec, —9.24°)  giving  sfn  =  0.8913,5^21  =  I  826 
and  All  =  —0.6709.  In  region  III,  the  Volterra  series  has  been  expanded  about  the  e  lui- 
librium  point  (ao,9o,<5eo)  =  (17.6°,  — 7.911239‘’/sec,  — 1 1.4°)  giving  gw  =  —0.180701,(721  = 
-0.361403  and  An  =  0.459482.  The  calculations  for  gw,g2\  and  An  for  region  II  are  giver 
in  Appendix  C.  In  the  prediction  analysis  we  use  the  following  procedure.  Starting  from 
conditions  q°  =  11°, 9°  =  0  and  equal  to  a  constant  value,  the  Volterra  model  of  the 
first  subspace  describes  the  airplane’s  response  until  a  =  14.36°  is  reached.  At  this  point 
the  Volterra  model  of  the  second  subspace  describes  the  response  until  q  reaches  one  of  the 
boundaries,  14.36°  or  15.6°.  If  the  response  intersects  the  14.36°  boundary,  then  the  first 
Volterra  model  governs  the  response.  If  the  response  intersects  the  15.5°  boundary,  then 
the  third  subspace  model  governs  the  motion  while  q  remains  between  15.6°  and  19.6°.  The 
fourth  subspace  model  governs  the  model  when  a  is  between  19.6°  and  28°.  Upon  entering 
a  new  region,  the  equilibrium  state,  xq,  is  changed  to  correspond  to  the  equilibrium  point 
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used  in  this  region.  The  initial  linear  state  a?!  is  set  equal  to  x  —  xq,  the  value  of  the  state 
at  the  boundary  crossing  minus  the  new  equilibrium  value.  The  initial  “nonlinear”  states, 
X2,  X3, ...,  are  set  equal  to  zero. 


We  define  Xj  as  the  state  resulting  from  a  Volterra  approximation  with  j  terms 

1 

Xj  =  xq + = 

t=l 

where  Xq  is  the  equilibrium  state. 

In  the  example  described  above  we  have  j  =  3 

•*'3(0  =  ^'o  +  •'TiiO  T  -^2(0  T  •'^3(0  (3.15) 

where  x,(f),i  =  1, 2, 3  satisfies  (3.9).  The  solution  to  the  nonlinear  differential  equations  (3. 1 ) 
and  (3.2)  is  denoted  by  x(t)  =  [a(t),g(f)]^.  We  use  a  fourth-order  Runge-Kutta  routine  to 
calculate  x(0i  thus  a{t). 

We  consider  two  step  inputs  6e  =  —9.2°  and  6e  =  —9.4°.  The  onset  of  the  limit  cycle 
occurs  between  these  two  inputs.  For  6e  =  —9.2°,  the  responses  a{t)  and  Q!3(<)  are  compared 
in  Figure  3-3.  Their  differences  are  presented  by  the  solid  line  in  Figure  3-6.  We  observe 
no  limit  cycle  for  this  input.  The  three-term  Volterra  response  is  an  excellent  match  to  the 
nonlinear  response.  The  maximum  difference  is  only  0.05  degrees  as  shown  in  Figure  3-6. 

For  =  —9.4°,  the  responses  a(<)  and  Q’3(t)  are  compared  in  Figure  3-4  with  their 
differences  presented  by  the  solid  line  in  Figure  3-7.  This  input  of  6e  =  —9.4°  generates 
a  limit  cycle.  Again  the  three- term  Volterra  model  response  is  an  excellent  match  with  a 
maximum  difference  of  0.01  degrees  over  a  time  interval  of  100  seconds  as  shown  in  Figure 
3-7. 

For  comparison,  we  consider  a  piecewise-linear  model,  choosing  the  following  regions  and 
lineariz('d  models  of  ('\(a').  Figure  3-2, 

C,(q)  =  (-0.07281587)0,  q  <  14.74"  (3.16a) 

G(o)  =  -1.07.3305924  +  0.0884 70922(q  -  14.74°),  14.74°  <  a  <  17.4°  (3.16b) 


(3.14) 


C,(a)  =  -0.8308956  +  0.03309905(q  -  17.4"),  17.4"  <  a  <  18.87" 


(3.16c) 


C,(q)  =  -0.7882234  -  0.016633734(a  -  18.87"),  18.87"  <  q  <  28"  (3.16d) 

The  responses  using  the  above  approximations  are  compared  with  the  nonlinear  responses, 
a{t),  in  Figure  3-5.  The  differences  are  presented  by  the  broken  line  in  Figure  3-6  for 
8e  =  —9.2",  and  in  Figure  3-7,  for  Sg  =  —9.4".  The  maximum  difference  in  Figure  3-6  is 
about  2.5  degrees,  and  in  Figure  3-7,  about  0.7  degrees.  These  differences  are  one  to  two 
orders  of  magnitude  more  than  those  obtained  using  the  Volterra  approximation. 
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rv  (DEGREE) 


REGION  I:  a  <  14.74'’ 

REGION  II;  14.74'’  <  a  <  17.4'’ 

REGION  III;  17.4"  <  a  <  18.87'’ 

REGION  IV:  18.87'’  <  a  <  28'’ 

FIGURE  3-2.  PIECEWISE-LINEAR  FIT  (SOLID  LINE)  FOR  ^(a) 
CURVE  (BROKEN  LINE). 


FIGURE  3-3.  THREE-TERM  VOLTERRA  RESPONSE  (BROKEN  LINE) 
PROVIDES  EXCELLENT  APPROXIMATION  FOR  NONLINEAR  SYSTEM 
RESPONSE  (SOLID  LINE)  FOR  8,  =  -9.2“,q°  =  lU  AND  0“  =  0. 


FIGURE  3-4.  THREE-TERM  VOI.TERRA  APPROXIMATION  (RROKEN 

LINE)  PREDICTS  NONLINEAR  SYSTEM  (SOLID  LINE)  LIMIT  CYCLE 
FOR  =  Jl"  AND  -r  (). 
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FIGURE  3-5.  PIECEWISE-LINEAR  RESPONSE  (BROKEN  LINE)  FAILS 
TO  MATCH  NONLINEAR  SYSTEM  RESPONSE  (SOLID  LINE)  FOR  S,  = 


=  lU  AND  r/"  =  0, 


Eli-ROR 

(DEGREE) 


TIME  (SECONDS) 


FIGURE  3-6.  COMPARISON  OF  THREE-TERM  VOLTERRA  AND  PIECEWISt 
LINEAR  ERRORS  FOR  <5,  =  -9.2^a“  =  IT  AND  7°  =  0  SOLID  LINE  REP¬ 
RESENTS  THE  ERROR  03(0  -  o(/),  BROKEN  LINE  REPRESENTS  THE 
ERROR  FOR  THE  PIECEWISE-LINEAR  APPROXIMATION. 
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FIGURE  3-7.  COMPARISION  OF  THREE-TERM  VOLTERRA  AND 
PIECEWISE-LINEAR  ERRORS  FOR  6^  =  -9.4°,  q°  =  11°,  AND  q°  =  0.  SOLID 
LINE  REPRESENTS  THE  ERROR  agCO  -  a{t),.  BROKEN  LINE  REPRE¬ 
SENTS  THE  ERROR  FOR  THE  PIECEWISE-LINEAR  APPROXIMATION. 
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3.2  Wing  Rock  Limit  Cvcle:  Example  2 

Unsteady  aerodynamic  effects  at  high  a  generate  a  wing  rock  limit  cycle  phenomenon  on 
aircraft  configurations  incorporating  slender  forebodies  (e.g.,  F-5  and  X-29),  Nguyen,  et  al 
[21,  20]  and  Ekandon,  et  al  [30].  We  consider  a  mathematical  model  of  an  experimental  wind- 
tunnel  wing  rock  model  developed  by  NASA  Langley  Research  Center,  Nguyen.  Whipple  and 
[Brandon  [24].  The  following  equations  (3.17)  and  (3.18)  hohl  for  a  vvmd-tunnel  sting  mounted 
model  on  an  apparatus  which  allows  the  model  to  rotate  freely  about  its  roll  axis  with  no 
angular  limitation. 

•p  =  p  (3.17) 

where  p  and  p  are  the  roll  angle  in  radians  and  roll  rate  in  radians  per  second,  respectively. 
Here,  the  constants  q,S,b,Ii  and  V’  are  the  dynamic  pressure,  wing  reference  area,  wing 
span,  roll  moment  of  inertia  and  free  stream  air  speed,  respectively.  The  coefficient  Cig  is 
the  rolling  moment  stability  derivative  due  to  sideslip  0.  The  coefficient  C/p  is  the  rolling 
moment  derivative  due  to  roll  rate  p  and  sideslip  rate  8.  Approximate  wind-tunnel  values  of 
these  coefficients  are  given  by  Nguyen,  Whipple  and  Brandon  for  the  angle  of  attack  a  =  30^. 


See  the  development  in  Appendix  D. 

Q^(q  =  ,30^)  = -0.4584  (3.19) 

C,p(Q  =  W\8)  =  0.4[1  -  7.64|/B|],  (3.20) 

0  <  1,^1  <  0.35 

p  =  28  (3.21) 

Eivaluating  the  constants  in  (3.18)  we  have,  see  [23], 

\  (^)  "  -26.6667  (3.22) 

(^)  "  0.764S5[1  -  3.82|^|]  (3.23) 
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In  view  of  (3.22)  and  (3.23),  Equation  (3.18)  becomes 

p  =  -26.6667(^  +  0.76485(1  -  3.82l<^l|p 

Our  nonlinear  model  is  described  by  (3.17)  and  (3.24). 

The  nonlinear  equation  (3.24)  lends  itself  to  separation  into  two  regions:  Region  1  with 
yP  >  0  and  Region  II  with  ^  <  0. 

p  =  -26.6667s?  +  0.76485(1  +  (-1  )^3.82v?]p  (3.25) 

where  j  denotes  the  region. 

Let  j"  =  (:^,  p]^  be  the  state  of  t  he  nonlinear  equations  (3.17)  and  (.3.24 ).  Let  x,  —  p,]  ^ 

be  the  contribution  to  the  state  from  the  ith  term  in  the  \'olterra  series.  The  differential 
equations  for  the  five-term  Volterra  series  approximation  have  the  form 


i]  =  /1ji  (3.26a) 

j2  =  -4X2  +  itl(3-l),  X2(0)  =  0  (3.26b) 

X3  = /1j3 -t-fif2(xi,X2),  X2(0)  =  0  (3.26c) 

x^  = +^3(xi,X2,  X3),  X4(0)  =  0  (3.26d) 

X5  =  Axs  +  g4(xi,X2,i3,X4),  X5(0)  =  0.  (3.26e) 


The  equilibrium  point  of  (3.17)  and  (3.24)  is  ipo  =  0,  po  =  0.  The  A  matrix  in  (3.26)  is  given 

4  _  [  0  1  1 

-26.6667  0.76485  ’ 

L  J 

The  g,  functions  in  (3.26)  have  the  form 

p.(x,,...,x,)  =  (-l)^2.92192p,2(x, . x.)  iV28i^) 

where  j  denotes  the  region  ( j  =  1  or  j  =  2)  and  where 

Pi2(xi)  =  .piP,  (3.28b) 
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522 (-^li  ■J'2)  ~  S^lP2  +  ‘P2P1 
532(-^’  1  ^3)  =  “plPa  +  <p2P2  +  >p3Pl 


(3.28c) 

(3.28d) 

(3.28e) 


542(Xi.  X2-  -r4)  =  'PlP4  +  ^2P3  +  ‘P3P2  +  <P4Pl 

.\s  in  cq.  (3.14),  wc  let  .r^  be  the  state  resulting  from  a  Vblterra  series  with  j  terms.  For 

j  =  1,2 . 5,  the  roll  angle  responses  for  the  initial  condition  <^(0)  =  5.73°  and  p(0)  =  0.0  are 

presented  in  Fignre.s  3  8  through  3-12  respectively.  1  he  roll-angle  response  of  the  linearized 
system  .ri  =  .4.ri  is  compared  with  the  nonlineai  response  in  Figure  3-8.  The  linearized 
system  diverges  to  a  roll  angle  of  180°  after  0  seconds;  it  does  not  predict  the  stable  limit 
cyc!^.  The  limit  ''ycle  is  |)rcdict(xl  by  the  two-t<‘rm  X’olterra  appro.ximation;  the  amplitude  is 
short  by  about  7  degrees.  The  three-term  \blterra  approximation  overpredicts  the  amplitude 
by  about  5  degrees.  The  four-term  and  five-term  \’olterra  approximations  provide  good 
predictions  for  both  the  amplitude  and  the  period  of  the  wing  rock  limit  cycle,  Table  3.1. 
The  five-term  Vblterra  approximation  predicts  the  amplitude  to  within  0.31  degrees.  The 
nonlinear  system  response  was  calculated  using  a  fourth-order  Runge-Kutta  routine. 

3.1.  Comparison  of  Predicted  Wing  Rock  Limit  Cycle  Characteristics 


System  v.  Roll  .Angle 


Vblterra 

.Approximation 

F’eriod 

(Seconds) 

Amplitude 

(Degrees) 

two- Terms 

1.21  on 

28.15 

three- Terms 

1 .2-257 

40.60 

four- Terms 

1.21.500 

34.46 

five- Terms 

1.218.57 

35.65 

Nonlinear 

1.2187 

35.34 

System 

We  observe  that  both  the  period  and  amplitude  values  for  the  Volterra  approximations 
given  in  Table  3.1  approach  the  corresponding  value  for  the  nonlinear  system.  Let  E,J  — 
2.  3.  4. 5  denote  the  error  for  the  amplitud('  as  predicted  by  the  i-term  Vblterra  approximation. 
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It  follows  that 


1^'. 


+1 1 


\E}\ 


<.4,  ;  =  2,3.4. 


riius  we  might  suspect  quadratic  convergence  for  the  sequence  We  make  this 

observa.tion  for  this  one  example  only  and  note  that  it  is  based  on  numerical  results  for  four 


terms. 
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TIM  10  (SIOCONDS) 

FIGURE  3-8.  LINEARIZED  SYSTEM  (SOLID  LINE)  DIVERGES,  FAIL¬ 
ING  TO  PREDICT  STABLE  WING  ROCK  OF  NONLINEAR  SYSTEM 
WITH  ^0'  ^  o.t;!".//’  0. 
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FIGURE  3-9.  TWO-TERM  VOLTERRA  APPROXIMATION  (BROKEN 
LINE)  PREDICTS  STABLE  WING  ROCK  OF  NONLINEAR  SYSTEM  (SOLID 
LINE)  FOR  =  0.73, =  0  WITH  7.2”  AMPLITUDE  ERROR. 


0  10  20 


TIME  (SECONDS) 

FIGURE  3-10.  THREE-TERM  VOLTERRA  APPROXIMATION  (BRO¬ 
KEN  LINES)  PREDICTS  STABLE  WING  ROCK  OF  NONLINEAR  SYS- 
I  E.\I  (SOLID  LINE)  FOR  v?"  =  Ti  TL/)”  =  0  WITH  5  2(i"  AMPLITUDE  ERROR. 


TIME  (SECONDS) 

FIGURE  3-11.  FOUR-TERM  VOLTERRA  APPROXIMATION  (BROKEN 
LINE)  PREDICTS  STABLE  WING  ROCK  OF  NONLINEAR  SYSTEM  (SOLID 
LINE)  FOR  =  5.73, 0  WITH  0.9^  AMPLITUDE  ERROR. 
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FIGURE  3-12.  FIVE-TERM  VOLTERRA  APPROXIMATION  (BROKEN 
LINE)  PREDICTS  STABLE  WING  ROCK  OF  NONLINEAR  SYSTEM  (SOLID 
LINE)  FOR  ifP  =  5.73,  p”  =  0  WITH  0.3U  AMPLITUDE  ERROR. 


4.  Analysis 

In  this  section  we  consider  several  approaches  to  analyzing  the  nonlinear  limit  cycles 
considered  in  Section  3.  First  we  suggest  an  approach  for  analyzing  the  wing  rock  linrit 
cycle  that  is  based  on  the  Volterra  submodels  developed  previously.  Following  this,  the 
usefulness  of  classical  approaches  to  the  analysis  of  the  wing  rock  and  longitudinal  limit 
cycles  is  explored. 

4.1  Volterra  Series  Analysis 

In  this  section  we  apply  the  Volterra  series  subspace  technique  to  the  wing  rock  example 
given  in  Section  3.2.  In  particular,  we  obtain  a  two-term  representation  for  the  solution  of 


the  system 


(^  =  p 


p  =  -26.67<p  +  0.76485(1  -  3.82|(p|]p 

For  convenience  we  let  a  =  26.67,  b  =  0.76485,  c  =  3.826  and  define  the  2x2  matrix  A  by 

^=(  ”  ll- 

—a  6 

The  system  (4.1)  can  now  be  represented  by  the  matrix  equation 

X  =  Ax  -h  g{x)  (4.2) 

where  x  =  and  the  nonlinear  vector  valued  function  g  is  defined  hy  g{x)  =  [0,  — c|(plp]^. 

It  is  to  be  noted  that  the  unique  equilibrium  point  for  (4.2)  is  Xq  =  0.  We  assume  that 
system  (4.2)  has  a  solution  of  the  form  x  =  kxi  +  k^X2  +  ...  where  k  is  a  constant  and 

1  i  =  li2,...  .  Substituting  x  into  (4.2)  yields 


kx\  -f-  k  X2**'  —  A{^kx\  -f-  k'^X2  -f-  ...)  -f- 


-c[kx\\  -|-  k^X2\  +  ...)(A:j:i2  +  k^X22  T  ••.] 


in  the  region  >  0.  At  this  point  we  equate  the  coefficients  of  the  first  and  second  powers 
of  the  constant  k  to  obtain  the  following  equations  for  the  first  and  second  terms  of  the 
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Vblterra  series  representation 


X]  =  Axx 


x-i  =  .4X2  +  S?i(xi) 


where  5fi(xi)  =  [0,-0x11X12]^.  To  include  initial  conditions,  say  x(0)  =  x°,  in  system  (4.2) 
we  modify  (4.3)  through  (4.4)  to  become 


Xi  =  >lxi,  Xi(0)  =  x*^ 


X2  =  Ax2 +5l(Xi),  X2(0)  =  0 

The  state  transition  matrix  for  (4.5)  is  given  by  (see  Appendix  F) 

^At  _  [  ^11  (^) 


Boiit)  B22(t 


where 


^ii(0  =  cosuot  —  B\2{t)  =  —sinuot 

2u>o  u;o 

—  a  n  /  \  ^  • 

“21(0  =  - SinujQt^  «22(0  =  COSLDot  +  - - SlUUJot 

ujq  2u;o 


and  uq  =  y  (4a  —  ^^)/4.  For  x°  =  [0,p°]^  we  have 


^Jf\-f,At  0  _  0  6(/2  Suit 

j  -  P  e 


X2[t)  =  e' 


0  +  Bx2{T)B22{T)dt. 


Using  various  integration  techniques,  we  find  that  X2(0  has  the  form  (see  Appendix  F) 


•28 


+  36a;o  —  24a  j  <1^52  3(^2^  ^ 

Combining  (4.7),  (4.8)  and  (4.10),  we  have  a  representation  for  the  approximate  solution 
X2{t)  =  xi(t)  +  X2(i)  for  (4.1)  with  initial  condition  (^(0)  =  0,  <^(0)  =  To  calculate  the 
amplitude  and  period  for  the  limit  cycle  of  the  wing  rock  example,  we  can  set 

X2(ti)  =  ari(<i)  +  X2(ti)  =  _p0  (4.11) 

where  tj  represents  half  the  limit  cycle  period.  Assuming  that  the  limit  cycle  is  symmetric 
with  respect  to  the  axis,  equation  (4.11)  allows  us  to  (at  least  in  theory)  solve  for  p°  and 
<1  in  terms  of  a,  b,  c. 

The  exact  representation  for  p°  and  in  terms  of  the  values  a,b  and  c  probably  can¬ 
not  be  obtained  in  closed  form,  and  ways  to  simplify  the  equations  to  obtain  approximate 
sol  .tions  should  be  considered.  From  this  example  we  see  that  simple  (two-term)  Volterra 
approximations  lead  to  complicated  expessions  when  solved  explicitly.  This  indicates  that 
techniques  of  a  more  qualitative  cind  indirect  nature  may  be  needed  to  show  how  the  explicit 
Volterra  solutions  can  be  simplified. 
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4.2  Classical  .Analysis 

In  this  section  we  present  several  cla.ssical  approaches  to  the  analysis  of  the  examples 
discussed  in  this  report.  A  review  of  several  standard  techniques  will  serve  to  put  the  new 
approach  we  have  considered  into  proper  perspective  and  to  suggest  methods  for  carrying 
out  analysis  using  Volterra  submodels.  In  the  following  sections,  the  wing  rock  and  longitu¬ 
dinal  limit  cycles  are  analyzed  using  an  energy  (Lyapunov)  method,  harmonic  balance,  and 
piecewise-linear  analysis. 

4.2.1  Wing  Rock 

The  equation  that  described  the  wing  rock  considered  in  this  report  is  of  the  form 

9  +  <np  —  Up  +  ^1^19  —  0  (4-12) 

where  p  is  the  roll  angle  in  radians  and  in  our  case  a  =  26.7,  b  =  0.765,  and  c  =  2.92. 

Most  classical  approaches  to  analysis  of  nonlinear  dynamic  behavior  assume  that  the 
behavior  is  characterized  by  a  single  nonlinear  differential  equation  (i.e.  one  region).  Since 
this  is  the  case  with  our  simplified  description  of  wing  rock,  it  is  an  ideal  candidate  for 
application  of  classical  methods. 

Energy  Analysis 

A  standard  equation  from  mechanics,  which  involves  nonlinear  damping,  is  of  the  form, 
see  [34], 

X  +  h{x,  i)  -{■  g{x)  =  0. 

The  energy  of  this  system  is  given  by 


and  the  change  in  energy  with  respect  to  time  by 

dE/dt  =  — x/if-T,  j). 

Our  e(|uation  (4-12)  can  be  put  in  this  form  using  the  following  correspondence 

X  =  p 


.40 


g  -  aif 


h  =  — -f 

To  evaluate  the  change  in  energy  of  our  system,  we  consider  two  cases. 

I.  >,?  >  0  :  — i)  =  (^^[6  —  C(^] 

=>  dEldt  >  (<)0  when  9  <  {>)bfc 

II.  (yp  <  0  :  — i/t(x,i)  =  (^^[6  +  C(^] 

=>  dE/dt  >  (<)0  when  if  >  (<)  —  b/c 

The  regions  of  energy  increase/dccrea.se  are  shown  in  Figure  4-1.  Ihe  situation  shown  in 
I'iguic  4-1  is  similar  to  that  of  the  classical  Van  der  Pol  oscillator.  I  he  origin  is  an  unstable 
eiiuilibrium  point  as  the  system  energy  increases  in  a  neighborhood  of  the  origin. 


FIGURE  4-1.  Regions  of  Energy  Increase/Decrease 


31 


As  the  system  energy  increases,  the  system  is  forced  into  the  region  \ip\  >  bjc  where  the 
energy  is  decreasing  and  in  turn  is  forced  back  into  the  region  \^p\  <  bjc  where  the  energy  is 
increasing.  This  alternation  between  increasing  and  decreasing  energy  results  in  the  observed 
limit  cycle.  If  the  coefficient  of  the  nonlinear  term,  c,  should  change  sign,  then  energy  will 
increase  everywhere  in  the  pha.se  plane  and  the  system  will  be  unstable. 

This  analysis  also  indicates  why  the  linear  submodel  approach  fails  to  predict  a  stable 
limit  cycle.  The  linear  submodels  accurately  redect  the  behavior  in  a  neighborhood  of 
the  origin,  where  the  energy  is  always  increasing,  and  thus  predict  completely  unstable 
behavior.  The  Vblterra  submodels  reflect  the  behavior  over  a  larger  region  of  the  state 
space  and  somehow  capture  the  positive  damping  that  occurs  for  larger  values  of  |<^|.  .Xii 
area  for  future  investigation  would  be  to  characterize  the  damping  in  a  truncated  Volterra 
approximation. 

Although  energy  analysis  is  a  very  powerful  tool  for  this  second  order  example,  its  use¬ 
fulness  declines  as  the  systems  become  more  complicated.  It  is  not  generally  possible  to 
define  a  useful  energy  (Lyapunov)  function  for  a  complicated  system  of  equations. 

Harmonic  Balance 

To  apply  the  first-order  harmonic  balance  method,  we  assume  that  the  oscillation  is  a 
pure  sinusoid 

V?  =  Acos{wt)  (4-13) 

of  unknown  amplitude  and  frequency.  Substituting  (4-13)  into  the  original  differential  equa¬ 
tion  yields 

—  Aw^ros{ivt)  4-  nAcos{ivt)  +  bAio.<>in(wl)  —  c/ln;lAco.s(?r/,)|.‘<(7;(?/’/)  =  0  (4-11) 

Substituting  the  Fourier  series  for  |cos(j/;<)| 

2  4 

|ros(?e<)|  =  — 1 - ro,s2?(d  +  .  .  . 

TT  37r 

into  (4-14)  and  dropping  all  harmonics  yields 

2c  .  .  4 

-  Air^rosi  u't)  +  a  Aros{  u'l)  +  bAwsi)i{trt.) - A^irsinl  ir! ) - rA^ivs77i(  —  irt)  =  0 

TT  fiTT 


Equating  the  coefficients  of  the  cos{wt)  tenns  gives  the  equation  for  the  limit  cycle  frequency 

w  =  \fa. 

Equating  the  coefficients  of  the  sin{wt)  terms  gives  the  equation  for  the  limit  cycle  amplitude 


Substituting  the  values  of  the  parameters  a,6,  c  for  our  case  gives 

te  =  5.17  (T=  1.216) 

A  =  35.37° 

which  compare  very  favorably  with  the  simulation  results  T  =  1.219  and  A  =  35.34°. 

The  harmonic  balance  method  is  a  very  powerful  tool  for  analyzing  simple  oscillations 
and  can  be  extended  to  higher  order  systems  with  a  concomitant  increase  in  complexity.  In 
the  next  section  we  apply  this  technique  to  a  more  complicated  example. 

4.2.2  Longitudinal  Limit  Cycle 

An  important  feature  of  the  equations  modeling  the  longitudinal  limit  cycle  is  that  the 
nonlinearities  cannot  be  represented  in  a  simple  manner  over  the  entire  range  of  interest. 
Most  classical  analysis  techniques  do  not  apply  in  this  situation.  The  approach  outlined 
in  this  report  decomposes  the  state  space  into  regions  in  which  the  nonlinearity  can  be 
represented  by  a  low  order  polynomial.  This  section  examines  two  other  approaches. 

The  first  approach  approximates  the  nonlinearity  by  a  piecewise-linear  function  and 
provides  qualitative  information  about  the  behav'or  of  the  resulting  system.  The  second 
approach  represents  the  nonlinearity  numerically,  using  a  spline  fit  to  wind  tunnel  data,  and 
investigates  the  usefulness  of  harmonic  balance  techniques. 

Piecewise  Linear  Analysis 

VVe  consider  the  longitudinal  equations  (3-1),  (3-2)  and  the  piecewise-linear  approxima¬ 
tion  of  Cz{a)  given  by  (3-16).  The  system  is  unstable  in  regions  II  and  III  (14.74°  <  o  < 
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18.87’),  as  the  coeificient  of  a  in  ("\(a)  is  positive,  and  stable  in  regions  I  and  IV.  For  a 
given  constant  input.  6^^,  the  equilibrium  value  of  o  is  found  by  setting  q  =  0  and  given  by 

«O  =  0.5-  1.5<5eo. 

For  values  of  between  0"  and  9.5°,  the  equilibrium  value  of  a  is  in  region  I.  Since  this 
region  is  stable  we  would  expect  a  stable  response.  F^or  values  of  S^q  between  —9.5°  and 
-12.2°,  the  equilibrium  value  of  a  is  in  an  unstable  region.  The  trajectory  cannot  converge 
to  this  equilibrium  point  and  it  is  forced  into  one  of  the  stable  regions  that  in  turn  forces  it 
back  towards  the  equilibrium  point.  Thus  we  would  expert  a  limit  cycle  to  develop. 

This  simple  analysis  allows  us  to  conclude  that  a  limit  cycle  will  develop  for  certain 
values  of  and  to  estimate  the  onset  of  the  limit  cycle  as  6^^  ^  —9.5°.  A  more  extensive 
development  of  this  analysis  can  be  found  in  [35].  In  addition,  it  should  be  noted  that  a  large 
amount  of  qualitative  and  quantitative  work  has  been  done  in  the  general  area  of  piecewise- 
linear  analysis.  One  advantage  of  using  Volterra  submodels  will  be  to  reduce  the  number  of 
regions,  which  may  make  analysis  more  tractable. 

Harmonic  Balance 

To  investigate  the  usefulness  of  harmonic  balance  techniques  in  more  complex  situations, 
we  consider  the  following,  more  accurate,  longitudinal  model  for  the  aircraft 

d  =  (7  +  0.1066  +  0.1097[C'^(a,6e)  +  C,,(q)7(0. 0123)] 
q  =  5.84[C„(a,,50  +  C,„,(a)7(0.0123)] 

where  a  is  in  radians  and  q  in  radians/second.  Because  of  its  small  magnitude,  the  ('’-^(a) 
term  will  be  ignored.  The  C\,Cm,Cmq  functions  will  be  rei)resented  by  spline  fits  to  wind 
tunnel  data.  Plots  of  these  functions  are  shown  in  Figures  3-1  (for  =  —18°),  4-2,  and  4-3. 


.3-1 


FIGURE  4-2.  T-2C  WIND  TUNNEL  DATA:  PITCHING  MOMENT  CO¬ 
EFFICIENT  FOR  ZERO  SIDESLIP. 
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FIGURE  4-3.  T-2C  PITCH  RATE  DERIVATIVE  DATA  FOR  ZERO 
SIDESLIP. 
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The  most  severe  nonlinearities  occur  near  a  =  15°,  which  is  the  onset  of  stall. 
Given  a  constant  8^,  we  first  assume  the  solution  has  the  form 


o(<)  =  7  +  acos{tvt) 
q[t)  =  7^  +  aqCos[wt)  bgSin(wt). 

Using  a  truncated  Fourier  series,  we  can  approximate  the  outputs  of  the  nonlinearities  by 

—  7m  T 

Se)  =  7z  +  a:Cos{wt) 

■  Tmg  “b 

where  the  7^  and  aj.  coefficients  are  computed  numerically.  Substituting  these  expressions 
into  the  model  and  retaining  only  first  harmonic  terms  results  in  three  equations  for  the 
unknowns  7,  a,  and  w.  The  harmonic  balance  (numerical)  solutions  for  a  and  w  are  compared 
with  values  extracted  from  time-domain  simulation  results  for  three  values  of  65  below 
Harmonic  First-Order  Balance  Simulation 


Se 

a 

w 

a 

w 

-15 

2.23 

4.71 

3.15 

4.16 

-20 

3.54 

3.56 

4.60 

3.19 

-25 

4.35 

2.87 

5.85 

2.73 

The  values  predicted  by  a  first-harmonic  analysis  are  not  very  accurate. 

One  problem  with  the  above  procedure  is  that  higher  harmonic  terms  in  the  Fourier 
expansions  of  Cmq{otii))  and  q{t),  which  were  ignored,  can  “beat  down”  and  affect  the  con¬ 
stant  and  first-harmonic  terms  through  the  Cmq{oi)^  term  in  the  q  equation.  By  considering 
the  multiplicative  interaction  of  the  second  harmonic  of  C'„,,(a(0)  with  the  first  harmonic 
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of  q  (only  one  of  many  higher-order  contributions)  the  improved  results  given  below  were 
obtained. 

Modified  Harmonic  First-Order  Balance  Simulation 


<^'e 

a 

w 

a 

w 

-15 

2.64 

4.46 

3.15 

4.16 

-20 

4.74 

3.21 

4.60 

3.19 

-25 

5.87 

2.76 

5.85 

2.73 

These  results  show  that  with  very  little  additional  numerical  complexity,  it  is  possible  to 
obtain  accurate  results  outside  the  stall  region.  Near  stall  ((5e  —15°),  however,  the  results 

degrade. 

A  more  accurate  (and  complicated)  analysis  assumes  the  solution  is  of  the  form 
oc{t]  =7  -f  aicos(tot)  -h  a2Cos(2u)t)  +  62sin(2u)f) 

-I-  a3Cos{3tvt)  -f  b^sinlSivt). 

.Assuming  similar  Fourier  expansions  for  q{t),  Cm.{oi{t),6e),  Cz(a(<),6e)  and  sub¬ 

stituting  these  expressions  into  the  model,  and  retaining  terms  through  the  third  harmonic 
results  in  seven  equations  for  the  unknowns  tu,  7,  Cj,  02,  ^3-  Solving  these  equations 

produced  the  following  results 

Harmonic  Third-Order  Balance  Simulation 


<5e 

a 

w 

a 

w 

-15 

2.96 

4.24 

3.15 

4.16 

-20 

4.36 

3.25 

4.60 

3.19 

25 

5.60 

2.7^ 

5.85 

2.73 

These  results  show  good  agreement  with  simulations  even  near  the  onset  of  stall.  From 
this  example,  accurate  predictions  of  the  amplitude  and  frequency  of  the  longitudinal  limit 


cycle  near  stall,  where  the  nonlinearities  are  the  greatest,  require  that  the  effect  of  higher 
harmonics  be  considered.  When  the  nonlinear  functions  involved  in  the  equations  are  derived 
from  wind  tunnel  data,  only  numerical  results  will  be  possible  in  general. 

In  summary,  the  classical  methods  provide  significant  insight  into  the  examples  consid¬ 
ered  in  this  report.  For  the  most  part,  however,  they  are  limited  to  assessing  stability  and 
predicting  limit  cycles,  and  therefore  do  not  form  a  broad  base  for  the  study  of  nonlinear 
flying  qualities.  On  the  other  hand,  the  Volterra  submodels  retain  the  essential  character¬ 
istics  of  the  nonlinear  model,  but  at  this  time  we  do  not  know  how  to  analytically  extract 
the  desired  information.  The  analytical  work  should  be  eased  somewhat  by  the  fact  that  the 
Volterra  approach  requires  fewer  regions  than  the  piecewise-linear  approach. 
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5,  Conclusion 

This  work  has  investigated  the  use  of  the  Volterra  series  representation  with  regard  to 
determining  nonlinear  flying  quality  parameters,  (also  see  [5,  36,  37]).  Our  first  goal  was  to 
determine  if  the  Volterra  series  could  accurately  represent  common  nonlinear  aerodynamic 
models  over  the  required  range  of  force/moment  conditions.  We  found  that  by  breaking 
up  the  state  space  into  natural  subspaces  in  which  the  nonlinearities  could  be  expressed  as 
low-order  polynomials,  the  total  nonlinear  model  could  be  represented  in  each  subspace  by 
a  low-order,  truncated  Volterra  series.  This  is  a  key  result,  since  dealing  with  the  original 
infinite  Volterra  series  would  be  intractable.  Our  simulations  showed  that  the  combination 
of  these  submodels  retained  the  essential  characteristics  of  the  underlying  nonlinear  model. 
Since  most  nonlinearities  encountered  in  practice  can  be  accurately  represented  locally  by 
p,olynomials,  this  technique  has  general  applicability. 

Our  second  goal  was  to  determine  if  a  piecewise-linear  model,  formed  by  replacing  the 
nonlinear  model  by  a  linear  model  in  each  subspace,  could  accurately  capture  the  nonlinear 
phenomena  contained  in  the  original  model.  We  found  that  this  was  not  the  case.  In 
simulations  of  a  wing  rock  limit  cycle,  a  piecewise-linear  model  did  not  predict  the  stable  limit 
cycle.  Using  two-term  Volterra  models  in  the  same  subspaces  did  predict,  via  simulations, 
the  existence  of  a  stable  limit  cycle.  In  simulations  of  a  longitudinal  limit  cycle,  the  piecewise- 
linear  model  was  significantly  less  accurate  than  the  low-order  Volterra  model. 

The  above  results  demonstrate  that  there  is  a  need  to  go  beyond  piecewise-linear  analysis 
of  nonlinear  systems.  One  appealing  feature  of  the  Volterra  series  is  that  the  “solution"  can 
be  written  down  explicitly,  as  shown  in  Section  4.1  for  the  wing  rock  example.  However, 
the  solution  contains  so  many  terms,  even  for  this  simple  example,  that  further  analysis  is 
necessary  in  order  to  extract  the  important  features. 

We  believe  the  first  priority  of  future  research  should  be  the  development  of  methods  for 
the  qualitative  and  quantitative  analysis  of  multi-region  phenomena  using  Volterra  submod¬ 
els.  If  such  methods  can  be  found  the  Volterra  approach  could  be  a  useful  tool  for  nonlinear 
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flying  qualities  research. 

Our  work  hcis  also  indicated  several  other  areas  that  would  benefit  from  further  research. 
The  first  area  is  the  selection  of  regions  and  the  trim  points  within  regions,  and  the  impact 
of  this  selection  on  the  accuracy  of  the  resulting  submodels.  In  the  examples  considered  in 
this  report,  the  selection  was  done  in  an  ad  hoc  manner.  The  second  area  is  the  connection 
of  submodels  at  the  region  boundaries.  There  are  several  ways  in  which  this  can  be  done 
and  the  method  used  in  our  examples,  initializing  the  first  (linear)  Volterra  term  with  the 
current  state  and  initializing  the  other  terms  to  zero,  may  not  be  optimum. 
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LIST  OF  SYMBOLS 


variables  units 


Description 


^e\ 


A 

b 

b 

Cm. 

c. 


9i-\ 

h, 

k 

P 

9 

9 

S 

u 


(deg) 

(rad) 

(deg) 

(deg) 

(deg) 

(rad) 


(ft) 

(1/rad) 

(1/deg) 

(1/deg) 


(rad/sec) 

(deg/sec) 

(Ibs/ft^) 

(ft^) 


Angle-ot-atta^'k 

Angle-of-sideslip 

Elevator  control  angle 

Equilibrium  elevator  control  angle 

Deviation  of  elevator  angle  from  equilibrium 

Euler  roll  angle 

Dynamics  of  linearized  equations  (i.e.  Jacobian) 

Input  matrix  for  control  variable  u  or  8e 
Wing  span  (Section  3.2) 

Rolling  moment  stability  derivative  due  to  sideslip  ^ 

Rolling  moment  derivatives  due  to  roll  rate  p  and  sideslip  rate  f) 
Pitching  moment  derivative  w.r.t.  a 
Pitching  moment  derivative  w.r.t.  6e 

Plunging  force  coefficient  (approximates  inverted  lift  coefficient) 
Nonlinear  term  of  ith  differential  Volterra  term 
ith  order  Volterra  kernel 

Grouping  term  parameter  in  Volterra  series  expansion 

Roll  rate 

Pitch  rate 

Dynamic  pressure 

Wing  reference  area 

Control  input 
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LIST  OF  SYMBOLS  (Continued) 


V  (ft /sec)  Free  stream  speed 

X  -  State  vector 

-  Initial  condition 

X.  -  State  vector  of  ith  term  in  Volterra  series  expansion 

X,  -  State  vector  of  Volterra  series  approximation  with  i  terms 
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Appendix  A 

Volterra  Series  Representation 

In  this  appendix,  we  show  how  the  Volterra  scries  arises  naturally  in  the  solution  of 
nonlinear  ordinary  differential  equations.  The  differential  form  of  the  series  that  we  have 
used  in  this  report  will  be  an  intermediate  step  in  the  calculation  of  the  integral  form  that 
has  been  used  by  Suchomel  [17]. 

To  keep  the  presentation  simple,  we  consider  the  scalar,  first-order,  nonlinear  differential 
equation 

X  =  X  bx^  -t-  u.  (A.l) 

The  procedure  we  outline  below  is  general,  and  can  be  extended  to  equations  of  any  order, 
see  Rugh  [20]. 

We  consider  the  solution  of  (A.l)  about  the  equilibrium  point  (xo,uo) 
considering  an  input  of  the  form 

It  =  A*«i, 

where  k  is  an  arbitrary  constant,  the  solution  can  be  written  in  the  form 

X  =  kxi  -h  k^X2  +  k^X'i  +  •  •  •  (A, 3) 

Substituting  (A. 2)  and  (A. 3)  into  (A.l),  and  equating  the  coefficients  of  the  various  powers 
of  k  (this  calculation  has  been  done  in  detail  for  the  specific  examples  considered  in  the 
report)  yields  the  set  of  equations 


j,  =  xi  -t-  U]  xi(0)  =  0 

(A.4) 

X2  =  X2  +  bx^  ^2(0)  =  0 

(A.5) 

=  (0,0).  By 

(A.2) 
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The  solution  of  (A. 4)  is  well  known  and  given  by 


3-1 


The  first-order  Volterra  kernel  is 


(A.6) 


hi{t  —  ai)  =  e‘  . 


By  considering  bxl  to  be  the  input  to  (A. 5),  the  solution  may  be  written  as 

Jo 

Substituting  for  Xi(s)  from  (A.6)  gives 

3^2(0  =  /  /  e*“‘''ui(<Ti)d(7i  /  e"~‘^'^u\{a2)d(T2ds. 

Jo  Jo  Jo 


To  facilitate  interchanging  the  integrals,  the  upper  limits  of  the  nested  integrals  will  be 
changed  from  s  to  t  by  inserting  an  appropriate  unit-step  function 


1  2  >  0 
0  2  <  0  ' 


The  result  is 
3-2(0 


=  f  *'6  f  j  <^_i(.S  —  '^'Wl(<3'i)dcr]  ,  6_i(s  —  I72)f  *  (CT2)dcr2d.s. 

Jo  Jo  Jo 

Changing  the  order  of  integration  yields 

3^2(0=  /  /  /  e*6_i(.s  -  crj)<5_i(.s  -  (72)d^!Ui(«7i  )u,((T2)d<T]d(72 

Jo  Jo  Jo 

The  functions  can  be  eliminated  by  writing 

rt  /•( 


.T2(0 


nbe'  f  c’dsti]{a\)U]{a2  .iaida2 
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Evaluating  the  innermost  integral  results  in 


-  ‘ ]  u  1  ( )  u  1  ( cr2 )  <71 


Thus,  the  second-oi  ler  Volterra  kernel  is 


To  put  this  in  the  time-invariant  form  h2(t  —  <Ti,<  —  <72)  we  can  write 


h2{t,(Tu(T2)  = 

=  h2it  -  <Ti,t  -  (72). 


Therefore,  the  solution  can  be  written  in  the  form 

X{t)  =  I  /li(<  -  (7i)u((7i)(f(7i  +  /  /  /l2(t  -  CTi,  t  -  (72)Ui((7,)Ui(<72)(f(7ic/cr2 

Jo  Jo  Jo 


which  is  identical  to  the  form  given  by  (2.1),  where  the  /?o(0  term  is  equal  to  zero  due  to 
our  selection  of  (0,0)  as  the  eciuilibrium  point. 
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Append i  ;  13. 

Longitudinal  Equations  of  Motion 

riio  follovvino;  longittidinal  equations  of  motion  in  wind  axes  foi  aircraft  flight  are  derived 
in  Ktkin  [2], 

0  =  q 

'  =  ~  -  o) 

'■*  =  9  +  -  /.]  +  ^cos(0  -  o) 

q  -  (^)  a  +  f^iLX  -  1,7.)  +  X[L:I\.  ~  l:i\) 
where: 

1.  The  states  O.\\o  and  q  denote  pitch  angle,  air  speed,  anglc-of-attack  and  pitch  rate, 
respectively. 

2.  The  quantities  m  and  I,  represen'  the  aircraft  mass  and  moment  of  inertia  about  the 
pitch  axis  y  through  the  center  of  gravity.  The  gravity  constant  is  denoted  by  g. 

3.  The  thrust  vector  in  wind  axes  is  denoted  by  {Tv,Ta)  where  Tv  is  the  thrust  along  the 
velocity  vector.  The  transformation  of  thrust  from  the  body  axes  representation  (7'j-.  T- ) 
to  the  wind  axes  yields  the  ecpiations 

I'v  =  TjCOf^a  +  7;,suna 

T„  =  —  7'j.sn7  0  +  T-cofia. 

4.  1  h('  aerodynamic  forces  in  wind  axes  are  denoted  by  drag  I)  and  lift  L.  These  are  related 
to  the  coi'fficients  f’/)  ^nd  f’/.  drag  and  lift,  respectively,  as  follows 


I)  =  Pioq  =  qSep 


L  =  Lift  =  qSCL 


where  q  is  dynamic  pressure  and  S  is  the  reference  area  of  the  wings. 

5.  Transforming  the  drag  and  lift  forces  from  wind  axes  to  body  axes  yields  the  relationships 

X  =  —Dcosa  -f  Lsina 

Z  —  —Dsina  —  Lcosa 

where  X  denotes  the  aerodynamic  force  along  the  body  j-axis  and  Z  denotes  the  aero¬ 
dynamics  force  along  the  body  2-axis. 

6.  The  aerodynamics  pitching  moment  coefficient  is  denoted  by  Cm- 

7.  The  mean  aerodynamic  chord  is  represented  by  c. 

8.  The  aerodynamic  forces  X  and  Z  and  the  thrust  forces  and  do  not  necessarily  pass 

through  the  center  of  gravity.  In  general,  they  have  the  moment  arms  and  /x,, 

respectively,  about  the  center  of  gravity. 

In  the  above  equation,  we  are  assuming  that  the  sideslip  angle  0  is  zero.  That  is.  we  are 
considering  pure  longitudinal  motion  in  the  vertical  plane  -  lateral  motion  is  considered  to 
be  zero. 

In  our  work  we  are  interested  in  the  cause  of  the  longitudinal  limit  cycle.  We  know  that 
it  is  generated  by  the  nonlinear  nature  of  the  lift  curve.  For  this  reason  we  simplify  the 
nonlinear  equations.  First,  we  consider  only  the  q  and  the  q  equations.  Second,  we  assume 
that  the  moment  arms  are  zero.  Third,  we  use  a  linear  model  for  the  pitching  moment  as  a 
function  of  only  angle-of-attack  a  and  elevator  control  Fourth,  we  use  a  linear  model  for 
the  elevator  control  6^  effect  on  lift.  According  to  the  assumptions  we  can  write 

Cl  =  Cifa)  +  Cu^^e 
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Cm.  —  Cmo  +  Cma<^  +  Cm^^^e 


Since  CL(a)  =  -  CDMtan{a) 


we  have  the  approximation 


Cl(q)  =  -C,(a)  +  C7,  10°  <  Q  <  20° 


where  C7  is  some  constant.  We  make  the  definitions 


Cl  = 


mV 


r  -  r 


To  ,  g  rn  ^ 


r 

M  —  T  ^rria 


c  -  ^C 
■*  1/ 


^ 

M3  —  »  ^mo  * 


The  resulting  equations  simplify  to  the  following: 


a  =  q  +  CiCz(a)  +  C2<5e  +  C3 


q  —  C4Q'  +  c^Sf.  +  Ce 

where  c, ./  =  1,2,...,()  are  treated  as  constants,  and  where  C\((y}  is  to  be  chosen  from  tin 
r('al  data  of  some  airplane.  For  the  T-2('  airplane  of  Ref.  22,  the  following  approximate 
model  is  obtained. 
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Q  =  <7  +  9.17r,,(Q)  -  1.8((5e  +  7°)  +  7.36 


9  =  5.73(C„„q  +  C„.,^6,)  +  2.9 

where  =  —1,  =  —1.5  and  where  the  nonlinearities  of  the  Cj(q)  curve  can  be 

represented  by  the  following  low-order  subspace  models: 


a  <  14.4°, 
14.4°  <  a  <  15.6°, 
15.6°  <  a  <  19.6°, 
19.6°  <  a  <  28°, 


C\(o)  ^  -0.0732a 
C,(a)  ^0.1  a^  -2.9a -h  20.0 
a(a)  ^  -0.02a2  -f  0.74a  -  7.8 
C,(q)  =  -0.47  -  0.02a. 


These  subspace  models  are  curve-fits  to  the  actual  post  flight  aerodynamic  data  of  the  T-2C 
obtained  by  system  identification  analysis,  [22] .  We  note  that  these  aerodynamic  derivatives 
represent  approximate  values  for  the  total  airplane.  The  above  model  may  lack  continuity 
at  the  boundary  points  since  all  numbers  have  been  rounded  off.  The  model  given  in  Section 
3  has  less  round-off  error  associated  with  it,  and  as  a  consequence  there  is  continuity  at  the 
boundaries.  This  continuity  is  with  respect  to  the  curve  but  not  necessarily  with  respect  to 
the  slope. 
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Volterra  Expansion  for  Longitudinal  Example 


In  this  appendix  we  present  the  details  of  the  Volterra  expansion  for  region  II  (14.36'’  < 
a  <  15.6°)  of  the  longitudinal  limit  cycle  example  considered  in  Section  3.1.  The  expansion 
for  region  III  is  similar,  and  the  equation  is  linear  in  the  other  regions. 


The  simplified  longitudinal  equations  are 

Q  =  q  +  9.168  a(a)  -  1.8336(^e  +  7°)  +  7.3619 
q  =  5.73(— Q  —  1.56e)  +  2.865. 

In  region  II  we  have 

C,(a)  =  0.09722o^  -  2.8653q  +  20.03846. 

Substituting  (C.3)  into  (C.l),  (C.2)  and  rearranging  gives 

•  \q  +  0.89130^  -  26.269Q  +  178.24]  ,  -1.8336]  r 

[  -5.73o  +  2.865  J  ^  -8.595 


(C.l) 

(C.2) 

(C.3) 

(C.4) 


where  x  = 

Considering  an  input  of  the  form 

6,  =  +  (C.5) 

the  solution  of  (C.4)  may  be  written  in  the  form 

.r  =  Xq  T  ^'-^i  T  k^X2  T  k^X:^  T  •  •  ■  •  (C.6) 

The  equilibrium  input  in  region  II  was  chosen  to  be  =  —9.24°.  This  resulted  in  the 

ecinilibrium  state  Xq  =  [1 4.36'’,  —  1 .7552°/.sec]^ .  Substituting  (C.5),  (C.6)  into  (C.4)  results 
in 
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kxi  +  k^X2  4-  k^X3  4-  •  •  ■  = 


(90  4"  kqi  4" .  •  •)  4"  .8913(q;o  4"  koti  4"  • .  •)^  —  26.269(q;o  4*  koi  4" 

— 5.73(q;o  4"  kcti  4" . .  •)  4"  2.865 


..)  4- 178.24 


,  —1.8336  fc  ,  If  \ 

-8.595  (®co  4- 


Using  the  equilibrium  values,  expanding  terms,  and  keeping  terms  up  to  order  3  in  k 


yields 


kxi  4-  k^X2  4-  k^X3  = 


A:(-0.6709ai  4-  -  1.8336^*, )  4-  A:2(-0.6709a2  4-  0.89130?  +  q^) 

lfc(-5.73a,  -  8.5956e,)  4-  P(-5.73o2) 


A:3(-0.6709a3  4-  1.7826aia2  4-  93) 

it3(-5.73a3) 


Equating  coefficients  of  equal  powers  of  k  yields  the  following  equations  for  the  Volterra 


scries  terms. 


-0.6709  1 
-5.73  0 


,  I  -1.8336  f 
^  '  -8.595 


-0.6709  1  ^  ,  0.8913a? 

-5.73  0  0 


-0.6709  1 


"  -5.73  0  0 


1.7826aiQ2 
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Appendix  D 

Wing  Rock  Model 

Aircraft  configurations  incorporating  slender  forebodies  generate  a  wing  rock  limit  cycle 
phenomenon  at  high  a  (angle-of-attack)  due  to  unsteady  aerodynamic  effects,  [23  -  33].  Delta 
wings  with  leading  edge  sweeps  greater  than  76“’  are  known  to  exhibit  wing  rock.  Nguyen, 
Whipple  and  Brandon  [24]  present  wind-tunnel  testing  results  of  such  an  80°  delta  wing. 
The  F-5  airplane  which  has  a  slender  forebody  is  known  to  exhibit  wing  rock  at  high  o. 
Lutze  [32,33]  presents  wind-tunnel  laterial-directional  aerodynamics  data  of  a  10-percent- 
scale  model  of  the  F-5E  airplane.  A  wind-tunnel  investigation  of  a  16-percent-scale  model 
of  the  X-29A  airplane  (which  has  the  Stmie  nose  section  as  the  F-5  airplane)  is  described  in 
Murri,  Nguyen  and  Grafton  [28].  Using  the  data  contained  in  these  references  we  construct 
a  wing  rock  model.  First,  we  investigate  the  aerodynamic  properties  of  the  80°  delta  wing 
that  is  inherent  in  wing  rock.  We  inspected  the  scaled  X-29A  wind-tunnel  model  and  found 
it  to  exhibit  similar  wing  rock  model  characteristics  as  the  80°  delta  wing.  For  this  reason 
we  construct  our  model  based  on  the  80°  delta  wing  data,  [24]. 

The  standard  wind-tunnel  test  technique  used  in  the  study  of  wing  rock  phenomenon 
is  the  free-to-roll  tests.  In  these  tests  the  physical  scaled  model  is  sting  mounted  on  an 
apparatus  which  allows  the  model  to  rotate  freely  about  its  roll  axis  (i.e.,  the  body  x-axis) 
with  no  angular  limitation.  A  80°  delta  wing  mounted  on  such  an  apparatus  is  shown  in 
Figure  2  of  Nguyen,  Whipple  and  Brandon  [24].  The  resulting  system  is  a  single  degree 
of  freedom  (SDOF)  system.  To  keep  our  presentation  as  simple  as  possible  without  loss  of 
qualitative  significance,  we  consider  a  mathematical  model  of  wing  rock  resulting  from  such 
a  wind-tunnel  testing  apparatus. 

In  a  free-to-roll  test,  the  pitch  angle  0,  which  is  fixed,  is  preset  according  to  the  desired 
angle-of-attack  (a)  at  zero  sideslip  [j3  =  0°).  The  roll  angle  ip  denotes  the  angle  of  roll  about 
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the  longitudinal  body  x-axis. 

Let  V  be  till'  tt)tal  idrsixx'd.  I'Ih'  eonventional  delinitioiis  of  tlu'  body  axes  velocity 
components  (u,v,w)  are 

u  =  Vcos{0) 
r  =  Vsin{0)sin{<f)) 
w  =  V  sin{0)cos{<j)) 


and  the  angle-of-attack  and  the  sideslip  angle  satisfy  the  identities 

,  .  10 
tanlQ)  =  — 
ii 

si7i{/3)  = 

From  these  we  have  the  identities 


Observe  that 


Vcos{0) 
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V 


cos{<f)\/V'^  — 


Since  u  and  \  '  are  constant  functions  of  time.  From  these  it  follows  that 


•  .  cos((^) + 

cos(/?)F 


or.  equivalently. 


,3  .  fo.s(<^)  . 

^  ^ 

cosyp) 


The  angle  of  attack  q  and  sideslip  angle  0  are  functions  of  0  and  if  and  satisfy  the  following 


relationships, 


tan(a)  =  tan(0)cos{f) 


(D.l) 


sin{3)  =  sin{9)sin{f) 


(D.2) 


tan{0)  =  sin{a)tmi{f) 


(D.3) 


cos{0)  =  cos{a)cos{0) 

.  ,  .  sin(0)cos{f) 

sin(a)  =  - — - 

cos{0) 

0  =  fsin(a). 


(D.4) 

(D.5) 


(D.6) 


In  our  construction  of  a  wing  rock  model  it  suffices  to  consider  the  following  approxima¬ 
tion  to  Equation  (D.3),  using  small  angles  in  0  and  f. 

0  =  fsin{a).  (D.7) 


The  roll  rate  p  satisfies 


f^p. 


(\).S) 


The  motion  about  the  roll  axis  is  governed  by  the  nondimensional  aerodynamic  rolling  mo- 


iiKMit  coefficicMit  (\  which  is  a  function  of  a,0,p,0  and  (S„: 

P -=  Cf{a,0,p.0,f^) 


(1).!)) 


.=>8 


where 


q  =  \p{h)V^  =  free  stream  dynamic  pressure 

p{h)  -  standard  air  density  at  altitude  h 

S  =  wing  reference  area 

b  =  wing  span 

li  =  roll  moment  of  inertia 

6a  =•  aileron  control  surface  deflection  angle. 

The  rolling  moment  coeflScient  Ci  has  the  expansion 

where 

0)  =  C'/p(a,  0)  +  Ci^[q,  /3)sin{a) 

and  where,  from  (D.6)  and  (D.8), 

Using  (D.7)  and  (D.IO)  we  rewrite  (D.9)  in  the  standard  form 

p  +  u;l{a)ip  +  2a»„(a)C(Q,v?)p  =  f{a,  l3)6a 


where 

‘^n(o)  =  -  Ctg{Q)sin{a) 

2u;n(a)C(o,<^)  =  -  =  ^sina) 


(D.IO) 

(D.lla) 

(D.llb) 

(D,12) 

(D.1.5) 

(D.14) 

(D.]5) 
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1)1  order  for  the  system  (D.8)  and  (0.12)  to  be  stable  at  high  a  it  is  necessary  that  the 
signs  of  C^,(o)  and  C(^[a,t3)  both  be  negative.  The  sign  of  Q^(q)  is  indeed  negative  at 
high  o.  We  remark  that  /^)  has  a  negative  sign  at  low  a.  But  airplanes  that  exhibit 

wing  rock  at  high  a,  25'’  <  o  <  50'’.  have  a  coefficient  derivative  C(^{a,l3)  that  is  positive  for 
small  sideslip  angles  and  negative  for  large  sideslip  angles,  Nguyen,  Whipple  and  Brandon 
[24]  and  Murri,  Nguyen  and  Grafton  [28].  The  wing  rock  limit  cycle  phenomenon  results 
from  this  switch  in  sign  of  Cfp(a,/^)  from  positive  to  negative  at  high  a  as  the  sideslip  angle 
increases  to  large  amplitudes. 


In  order  to  make  use  of  Eq.  (D.IO)  we  need  (quantitative)  expressions  for  the  aero¬ 
dynamic  derivatives  and  Wind  Tunnel  model  data  provides  C^j(a)  and 

Op(o.  d)  in  graphic  form.  We  remark  that  wind  tunnel  model  data  provides  a  model  for  the 
left-hand  side  of  Eq.  (D.Ua)  but  it  does  not  provide  separate  models  of  the  terms  fVp(o.d) 
and  CV,(a.d)  on  the  right-hand  side.  The  model  given  below  in  Eqs.  (D.16)  and  (D.17)  is 
an  approximation  to  the  wind  tunnel  data  plots  contained  in  [24|. 

The  model  of  Ctp{a,  d)  presented  in  Nguyen,  Whipple  and  Brandon,  [24],  for  the  80'’  delta 
wing  is  typical  of  airplanes  exhibiting  wing  rock  at  high  a.  Their  model  can  be  approximated 
by  the  following: 

NASA  WIND-TUNNEL  STING  MOUNTED  MODEL  80"  delta  wing 


0^0  =  30Nd)  =  0.4[I  -  7.64|d(],  0  <  |d|  <  0.35  (D.lfi) 

=  30'’)  = -0.4.584  (D.17) 

V  =  2d  (d  =  v^.suuo  at  a  =  30")  (D.IS) 

(~)  (d^)  ^  ^  =  .TOlSh]!  -  3.821|^|]  (1).20) 

(iO 


so  that 


P=  -26.6667V?  +  0.76485 [1  -  3.82|v?l]p  (D.21) 

riuMc  are  no  control  surfaces  on  their  wind-tunnel  model.  Therefore,  the  /{a,  (3)  term  drops 
out.  The  above  model  is  used  in  the  Volterra  Series  analysis  of  Section  3.2. 
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Appendix  E 

Expansion  of  Volterra  Series  for  Wing  Rock  Example 


The  model  for  wing  rock  as  derived  in  Appendix  D  is  given  by 


p  =  —26.7,^  +  0.76[1  —  3.82li^[|p. 


We  can  represent  this  equation  with  respect  to  the  following  two  regions: 


Region  I:  >  0 


p  =  —'26.1p  +  0.76li  —  3.S2>^]p 


Region  It:  o  <  0 


p  =  —26.7^  +  0.7G[1  +  3.82i^]p. 


We  note  that  the  p  equation  is  bilinear  in  the  ^pp  term. 


.A  state  space  representation  of  the  nonlinear  equation  is  given  by 


'■P  —  0  1  -L  ^  1  /  —  1  2 

p\  ~  [-26.7  0.76J  [p\  ^  [(-l)'2.92^pj  ’  '  “ 


Let  :r  =  ^  and  express  the  system  as 


.7  =  .-lx  +  .  ,  ,  .  c  =  ( - 1  )'2.92,  »  ^  1 , 2. 

(c)A /.  term  \ 


Using  the  Volterra  Series  technique  described  in  Section  2,  we  represent  the  nonlinear 
solutions  (pit)  and  p{t)  as  the  infinite  expressions 


p>{t)  —  p>o  P  +  •  ■  •  +  • 
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p{t)  =  Po  +  +  P2{t)k^  +  •  •  ■  +  Pn{t)k''  +  •  .  . 

(v’O'Po)  =  (0>0)  equilibrium  point  (unstable). 

Taking  the  product  of  >p{t)  and  p{t),  collecting  terms  with  respect  to  like  powers  of  k  and 
taking  into  account  the  equilibrium  point  yield 

ipp  =  {<-P\P\)k^  +  {^iP2  +  'p2P\)k^  +  {^\P3  +  ‘P2P2  +  ^3P\)k‘*  +  ■  ■  ■ 

We  make  the  identifications; 

Oui^i)  =  'PiPi 
P22(j'i,-'r2)  =  ‘P1P2  +92P1 
<733(^1 ,  a-2,  ^3)  =  <r>iP3  +  92P2  +  93P1 

i^P  —  9l2i^l)k^  +  922(^1^  ^2)k^  +  532(^1,  ^2>^3)k'*  +  .  • .  . 

The  Volterra  Se  nes  approximation  of  the  nonlinear  equation  is  therefore  given  by; 


xi  =  Ax\ 


X2  =  Ax2  +  Pi(a-i),  9iixi)  = 


0 

CPl2(.Tl) 


,  c  =  (-l)^2.92.  ;  =  1,2 


is  =  A.T3  +  (jr2(^i,^2),  92(^11^2)  = 


0 

(^922(^1^  ^2} 


in  —  -^^n  T  9n—l  (^1?  •  •  •  i^n— l))  9n—  \  (^Ij-'-i^n-l)  — 


0 

eP(n  —  1  )2 ( ^1  ^  •  •  ■  1  —  1 ) 


where  .r, 


V?. 

P. 


?'  = 
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Two-Term  Volterra  Series  Approximation 


In  tins  appendix  we  present  the  details  of  the  Volterra  series  analysis  for  the  wing  rock 
example  discussed  in  Section  1.1.  As  requested  by  the  contract  monitor,  all  derivations 
needed  to  obtain  the  two-term  Volterra  series  approximation  X2{t)  =  Xj(t)  +  X2(f}  for  the 
system  (4.1),  thus  (4.2),  are  presented  here  for  the  region  y?  >  0.  The  2x2  matiix  A  and 
the  constants  a.h,c  as  well  as  the  2  x  1  vector  functions  x,g  are  as  defined  in  Section  4.1. 
From  (4.5)  and  (4.6)  we  see  that  the  2  x  1  vector  functions  J’i(/)  and  X2{t)  satisfs' 


where 


.i'l  =  A,ri(/),  .r,(0)  = 


•''2(0  —  '‘Tr2(^ )  + 


0 

-r.ri,(/).ri2(/) 


.r2(0)  =  0 


Xi{t)  = 


a-nlO 


(F.l) 

(F.2) 


Solving  for  xAt)  in  (F.l) 


Fcpiation  (F.l)  is  a  first-order  linear  ordinary  differential  equation  with  constant  co¬ 
efficient  matrix  .4.  Figenvalues  for  the  matrix  A  are  Aj  =  {b  -|-  y/b'^  —  4a)/2  and  A2  = 
(6  —  yjb'^  —  Aa]/2.  Note  that  b^  —  4o  <  0  thus  Ai  =  (b/2)  +  iu>Q  and  A2  =  (6/2)  —  where 

a,'o  =  ( \/4n  —  6^)/2.  The  matrix  exponential  the  state  transition  matrix  for  (F.l),  can 
be  computed  by  the  identity  (see  Miller  [38,  page  111]) 

=  W\{t)l  -F  iV2{t){A  —  A] 7)  ( F.3) 

where  /  is  the  2x2  identity  matrix  and  the  scalar  functions  t/’i  and  W2  satisfy 

(■(',(/)  =  A|U>,(/),  "’i(O)  =  1  (FI) 

and 

=  Xzu'iil)  +  nq(0,  »'2(0)  =  0.  (F.5) 
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It  now  follows  that 


(F.6) 


and 

W2{t)  =  /  e~^^’’w2{s)ds 

Jo 

=  e^='‘  / 

Jo 

p^2t 

=  -i-p [£■<>■ -'>1  -  1),  (A,  -  A,)  /  0 

Ai  —  A2 

gAit  _ 

Ai  —  A2 

Using  identities  (F.6)  and  (F.7)  in  (F.3)  yield 


(F.7) 
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^  ^6t/2 


COSUlot 


bsiuLUot 


sinuot 


2u1q 
asinuot 

U>o 


UJo 


COSLUot  + 


bsinuot 

2a>o 


Let  denote  the  row-j^  column  entry  in  the  above  matrix,  see  equation  (4. 

The  solution  2:i(<)  for  the  initial  value  problem  (F.l)  is  given  by 


xi(0 


.At 


'  0  ■ 

'B,2{t)' 

p\ 

—  t  '  p 

B22\t) 

(F.S) 


Solving  for  .TtITI  in  (F.2) 


Equation  (F.2)  is  a  first-order  linear 
function 


ordinary  differential  equation  with  a  nonzero  forcing 


0 

-ca:i,(<)xi2(0 


where  from  above  xn(f)  =  and  Xi2(<)  =  e'’‘/V-^22(0-  The  variation  of 

constants  formula  together  with  the  zero  initial  conditions  of  (F.2)  yields 


X2(0  =  (P°)^  /  e"*'' 
Jo 


0 

— c 


e'^B,2(r)i?22(r)dr. 


We  now  employ  the  above  calculation  of  (in  terms  of  the  functions  defined  above) 

to  express  X2(0  by 

X2(0  =  (p“)"e^‘  /  [?1  (-c)e'’^B,2(r)522(r)dr 


0\7 „At 


rAf)  -  -c{pye 


'  B,2{-t)' 

L 

/?22(-r) 

Bi2iT)B22iT)dT 


6G 


X2{t)  =  5,2j-rj  B,2ir)B22{T)dT.  (F.9) 


The  matrix  is  nonsingular  and  has  an  inverse  e  We  shall  now  note  that  ^  0 
and  use  (F.9)  to  obtain 

^  j\6r/2  Bu{r)  B22{T)dT .  (F.IO) 

We  now  direct  our  attention  to  the  first  row  of  the  right  hand  side  of  (F.IO). 


Rn{t)  =  f  e>>^f^Bu{-T)Bi2{T)B22{r)dT 
Jo 

1  /"*  6  /"* 

= - 2  /  e^’^^^sin^UoTcosLOoTdT  —  —3  / 

^0  Jo  ^^0  Jo 


Integration  by  parts:  u  =  =  sinJujQTCosuiQTdr 

du  =  u  = 

2  ’  3wo 


—  — r  [  e^'^^^sinJu>oTdT 
2u;3  Jo 


Rii{t)  —  ~ 


e^’^^^sin^uot  b  /'  j^/2  3 

-■-  ...  —  -  f  4>  /  c»r> 

Jo 


^ sin  u2qt dr . 


For  the  second  row  of  the  function  on  the  right-hand  side  of  (F.IO)  we  have 


li2x(t)  =  f  e^^l^B22(-T)B22{T)Bn{r)di 
Jo 
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i?2i(0 

H2i(0 

^2l(0 


i: 


hrl'li  2  ^  ■  2  .SITIU^qT 

e  ^  {cos  luqT  —  - — -Sin  ujot) - dr 

^(i)^  (jJ(y 


I 


t  ^br/2  ^2 

- [1  —  (1  +  --^)sin^u2QT]sinLOoTdT 

0  ^*^0  4u.’q 


I 


t  gbr/2  M  g6r/2  ^2 

sinuJordr  —  I  - (1  +  — ■^)sin^UQTdT 

Jo 


Uq 


U)o  4Lt.’o 


We  now  integrate  the  first  term  on  the  right  to  obtain 


/?2l(f)  —  / .  2\ ~  UjQCOSuJQt)  +  u;o] 

UQ[b^  +  4a>5)  2 


rt  ^bT/2  p 

-  /  - (1  +  —-x)si'nJujQTdT. 

Jo  Wo  4u;^ 


In  order  to  eliminate  the  integral  terms  from  both  i'?n(0  and  R2\{t)  we  note  that  (see 
integration  tables) 

rt  4  6 

/  e^^^’^sin^Uordr  =  77— r^{(-smu;ot  +  -Zu)oCosuot)t^^/'^sin^ujQt 
Jo  0  +  dOWg  2 

+6wq  f  e^^l^sinuioTdT] . 

Jo 


Integrating  the  last  term  gives 


/ 


t^'^^^sin^iUoTdT  =  t:: - {{-sinuJot  ~  2u)oCosuJoi)^^^^^sin}ijJot 


62  +  36u;2  1^2 


24wo  u,/2/  ^  ■ 

+  .  .--W  ^  (-woro.swo/.  +  -.sjnwoO 


62  +  Iw,^ 
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(Note:  ^  =  4rt) 

4 


24u;g 

a{b''  +  36u;o) 


0(6^  +  36a;Q) 


We  now  have  the  following  representations  for  R\\{t)  and  /?2i(0- 


^hin 

Rii{t)  =  --:^sin^> 
SWq 


.  2  .  ,  24uJaCosuot  12uj^bsinu;ot .  24u;n, 

I2uostn^ujotcosu:ot - ^ ^  + - - - 

a  a  a 


■^21(0  =  - [e^^^^{-sinujot  —  u)ocosu)ot)  +  Wo] 

Qu-’o  2 


-  k  (sqrW) 


—  \2u}Qsin^ijjQtcosu}Qt  — 


24u}QCosu>oi 


12u!nbsinu}ot  ^  24a;n. 

+ - ^ - -)  + - -}. 

a  a 


It  follows  from  the  above  representations  for  Rii{t),  /?2i(0i  (F  IO)  that 


x,(t)  =  -e'"(p”)=c 
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-c(p°)V‘{ 


-(362  ^  3g^2) 
—Sab 


sin^ujot 


\  3ujo{b‘^  -f  SGwq) 


4b  sin^uotcosLOot 

12«J 
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-62  -  36a;2  +  24a 


COSLu'oi  ^2 

(6^  +  36uJo)a 


-862 

6^  +  366a;o2  -  24a6 


sinujQt 

2aa;o(62  +  36a>2) 


ef>(/2 


Appendix  G 

SMP  Function  for  Volterra  Equations  Calculations 

The  SMP  function  used  to  generate  the  differential  equations  for  the  Volterra  series  terms 
is  listed  in  Figure  Al.  The  nonlinear  differential  equation  is  assumed  to  be  of  the  form 

i  = 


riie  inputs  to  the  function  are; 

S/  =  {/i,/2- ••,/«}  =  components  of /(x,u). 

$var  =  {.Cl, . . . ,  u}  =  list  of  state  variables  and  the  input  variable  (last). 

$eq  =  list  of  equilibrium  values  of  the  variables  in  $var. 

%ord  =  number  of  Volterra  series  terms  to  be  generated. 

A  sample  output  of  the  program  for  the  wing  rock  example  is  shown  in  Figure  A2 

This  function  implements  the  procedure  outlined  in  Section  II.  A  line-by-line  description 
of  the  function  follows. 


1.  %n  is  set  equal  to  the  order  of  the  system. 

2.  %sub  is  set  equal  to  the  substitution  list. 

%ord 

%vari  —*  $uar,(jf](%a)-'  +  $eqi  f  =  1, . . . ,  %n 

j=i 

$uar%„+j  ->  (%a)$uar%„+,  +  $eg%„+i 


3.  %fs  is  a  set  equal  to  $/  with  the  above  substitutions. 

4.  Technical  detail.  a  projection  that  converts  truncated  series  expansions  into  poly¬ 

nomials,  is  forced  to  distribute  over  lists. 
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5.  The  function  %fs  is  expanded  in  a  Taylor  series  with  respect  to  the  dummy  variable  (%a) 
up  to  order  %ord. 


b.  For  each  V’olterra  series  term  is  set  equal  to  the  coefficient  of  (%n)'  in  the 

expansion  of  % f  s. 


T1 


r 


/*ivoltt$f,$var,$eq,$ord] 

gives  differential  equations  for  the  components  of  the  state 
corresponding  to  terms  in  the  volterra  series  for  the  state 
equation  dx/dt=$f (x,u) ,  The  state  variables  and  input  variable 
are  listed  in  $var  with  the  input  variable  listed  last.  $eq 
is  a  list  of  the  equilibrium  values  of  $var,  and  $ord  is  the 
number  of  terms  in  the  volterra  series  to  be  considered.  M/ 

volt($f,$var,$eq,$ord]i i\ 

( ';n  :  L  en  t  $  f  1 ;  V 

:';sub:Cat[Ar[Xn.$vartSl]->Sum[[$var[$l]][i]  '/.a^i  ,  {i,l,$ord}]S 
+  $eq[$l]],  {$var[%n+l  ]->/Ca  Svart?in+1  ]  +  $eqt5Cn+l  ]>  ]  ;S 
%fs!St$f.Xsubl;_Ax[Ldist]:l;5lfe:AxIPs[/ifs»?Ca,0,$ord]];\ 
Do[i.l.$ord,Dotj,l,Xn,Pr["d>'dt”,[$var[j])ii],”  =  ",\ 
N[Col[Ex[Coef[:ia'i,:{fe[j]]ni]];Prn;Prt)]) 


FIGURE  A*.  SMP  Function 

SMP  1.5.0 

23-JUL-1987  08:57iA1.05 

•  HI]::  f :  {p,-26.6667  phi+. 76485  (1-3.82  phi  )  p) 

»0[1]:  {p, -26 . 6667phi  +  0.76485p  (1  -  3.82phi)} 

11(2]::  vol  terra 


voltCf, {phi,p,u}, (0,0.0}, 51 


d 'dt  phi 1 1  ] 

d/dt  pCn 


d/dt  phi [ 2  1 

d/dt  p[2] 


d/dt  phi[31 

d/dt  p [ 3 1 


d/dt  phi (41 

d/dt  p(4] 


d/dt  phi [51 

d/dt  p[5] 


ptll 

0.76485p(ll  -  26 .6667phi(l ] 


p(21 

0.76485p(2]  -  26 .6667phi[2]  -  2.92175ptl]  phidl 


p(3] 

0.76485pr33  -  26 . 6667phi [ 3 1  -  2.92173p[n  phi[2] 
-  2.92173pt2]  phitll 


p[4] 

0.76485p(43  -  26.6667phil4]  -  2.92173p[l]  phir3] 
-  2.92173p(2J  phi[2]  -  2.92173p[3]  phi [ 1 1 


pC5I 

0.76485p(53  -  26 . 6667phi  [  53  -  2.92173p[13  phi[43 

-  2.92173p(23  phi  ( 3’3  -  2.92173p[31  phi[23  . 

-  2.92I73p[43  philll 


#1(43::  Exitci 


FIGURE  A2.  OUTPUT  FOR  THE  WING  ROCK  EXAMPLE 
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